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We introduce a computational framework which avoids solving explicitly hydrodynamic equations 
and is suitable to study the pre-merger evolution of black hole - neutron star binary systems. The 
essence of the method consists of constructing a neutron star model with a black hole companion 
and freezing the internal degrees of freedom of the neutron star during the course of the evolution of 
the space-time geometry. We present the main ingredients of the framework, from the formulation 
of the problem to the appropriate computational techniques to study these binary systems. In 
addition, we present numerical results of the construction of initial data sets and evolutions that 
demonstrate the feasibility of this approach. 
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I. INTRODUCTION 

The coalescence of compact binaries are among the most important sources of gravitational waves to be detected 
by ground-based laser interferometers like LIGO, TAMA, GEO, and VIRGO. These observational efforts have served 
as the primary motivation to simulate compact binary systems. The simulations will provide crucial information in 
support of the data analysis. Most of the focus has been on the simulation of black hole (BH) binaries and neutron 
star (NS) binaries, leaving aside the BH-NS binary system (see 1, 2, 3] for recent progress). From the point of view 
of Numerical Relativity, the system BH-NS exhibits a dual challenge. It has the difficulties of evolving the geometry 
in the vicinity of BHs together with the difficulties inherent to magneto-hydrodynamical calculations. 

From an astrophysical point of view, NS-BH binaries have the added interest of their potential relevance to gamma 
ray bursts. Observations of short gamma ray bursts (see 0, IE IE an d more recently [EIEEJ) suggest that the 
coalescence of NS binaries and/or BH-NS binaries are the underlying mechanism in the central engine for those bursts 
that last about 0.3 s . 

Due to the high computational cost of simulations of BH-NS binary systems (and the high cost of developing the 
appropriate computational infrastructure) it is desirable to have a framework in which to study certain dynamical 
regimes of this system, relevant to gravitational- wave physics, that admit the use of approximations and thus facilitates 
a drastic reduction of the computational resources required. In this sense, we are particularly interested in the stages 
of the evolution where the predictions of post-Newtonian methods become unreliable and where numerical relativistic 
simulations should take over. The hydro-without-hydro approximation that we propose in this paper aims at covering 
such part of the dynamical regime where the dynamical timescales related with deformations of the NS due to tidal 
effects are much bigger than the orbital timescales. In this regime we can freeze most of the hydrodynamical degrees 
of freedom and evolve a finite number of them by using appropriate approximations, focusing the attention on the 
radiation reaction effects in the orbit. It is this reduction of degrees of freedom which avoids the use of hydrodynamical 
computations. Such an approach may also be relevant for extreme-mass-ratio binary systems whose dynamics are 
driven by radiation-reaction effects. The main goal of this paper is then to describe the setup for this type of 
simulations and test, by extending an existing code for BH evolutions, the basic ingredients of this framework. These 
ingredients can be summarized as follows: 

a. Spacetime description: We use a full numerical relativity description of the spacetime based on the 3+1 BSSN 
formulation (Baumgarte, Shapiro [TTj . and Shibata, Nakamura |12 |) of Einstein's equations as implemented in the 
3D numerical relativity code MAYA (described in detail in |l3i Il4j). The MAYA code is written by using the 
computational toolkit CACTUS [r| and the Mesh Refinement package CARPET [IE 113 • The key extension of this 
code regarding this project is the introduction of the matter source terms associated with the NS and the associated 
infrastructure to evolve the matter degrees of freedom. 
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b. Construction of Initial Data: Initial data is constructed by su perpos ing; a Schwarzschild BH in Kerr-Schild 
(KS) coordinates with a boosted Tolman-Oppenheimer-Volkoff (TOV) model in a coordinate system that has 
similarities with the ingoing-Eddington-Finkelstein coordinate system. The TOV model describes the equilibrium 
configuration of a stationary, relativistic star in spherical symmetry and thus represents a convenient description of 
the matter sources in the BH-NS system at hand. The superposition is carried out by using a procedure analogous to 
that used for the construction of binary BH initial data initial data in KS coordinates 20] . Even though the resulting 
initial data do not satisfy the constraint equations exactly, the constraint violations are comparatively small because 
the individual BH and NS models do represent exact solutions of the Einstein field equations. In this work we will 
investigate in detail the properties of these initial data from a physical and mathematical point of view. We will 
also discuss potential modifications of the data with regard to their repercussions on the physical content. Finally 
we illustrate how adequate evolutions, even in the case of more extreme mass ratios, are facilitated by the modern 
infrastructure of numerical relativity. 

c. Description and motion of the NS: In the hydro without hydro approximation pursued in this work, we need 
to address two issues relating to the description of a particular type of energy-momentum distribution: the update 
of the matter sources that appear in the BSSN evolution equations and the motion of such a distribution. In a full 
numerical relativity setup this is not necessary as the complete evolution of the matter sources is taken care of by 
the Einstein field- and the energy-momentum conservation equations. In reducing the degrees of freedom in the way 
we are proposing, one is usually left with an approximation scheme in which the internal and external motions are 
cleanly separated, even though they are coupled in general. The external motion consists in the motion of a reference 
point of the matter distribution, typically a relativistic generalization of the Newtonian concept of center of mass, 
whereas the internal motion consists in the evolution of the parameters describing, for instance, the deformations of 
the matter distributions due to tidal deformations. As a consequence of reducing the number of degrees of freedom, the 
description of the motion of the energy-momentum distribution in our setup consists of a set of ordinary differential 
equations (ODEs). The matter sources, as they appear in the evolution equations for the spacetime geometry are thus 
updated by merely substituting the values of the energy- momentum quantities after the motion of the matter has been 
determined. In this paper, we consider the simplest such description of the NS, namely that of a star rigidly moving 
along a prescribed trajectory. By using this simplification we ignore internal motions as well as radiation-reaction and 
finite size effects due to the external gravitational field. More specifically, the matter sources are prescribed by taking 
particular density and pressure profiles from a TOV model which are going to be maintained along the evolution. The 
only way in which the sources change is through the rigid bulk motion of the center of the TOV model along a fixed 
trajectory. The reason for considering such a simplified model is to perform a feasibility study of this method and 
thus lay the foundation for future studies of more realistic physical configurations of the BH-NS binary system. A 
similar approach to that proposed here has recently been used by Bishop et al pi) in the framework of a characteristic 
formulation of the Einstein equations. 

The plan of this paper is as follows: In section|H]we present the mathematical description of the different ingredients 
of our framework: spacetime geometry and dynamics (|II Afl . the descriptions of the BH and the NS l|II Bl and III CI 
respectively), the initial data for the BH-NS system Ijll Dfl . the dynamics of the NS (|II E|l . In section llll Al we describe 
the computational framework and report on the results of numerical tests we have performed to analyze the initial 
data (subsection IIII B|l and to test time evolutions with matter sources l|III C|) . We conclude in section llVl with a 
discussion of our findings. For the purpose of presenting the relativistic equations we set c = G = 1 throughout this 
work. 



II. THEORETICAL FOUNDATIONS 

In this section we describe the different ingredients that constitute the mathematical framework that we are propos- 
ing for the numerical description in the time domain of BH-NS binary systems. For this purpose we need to consider 
both, the generation of initial data, and the evolution of these data. In the construction of initial data we closely follow 
a procedure used for a binary BH systems: the superposition of two BH solutions in KS form (cf. |20l l22l 123. l24|). 
This technique is based on the properties of the KS form of a single BH metric and its invariance under Lorentz 
transformations of the coordinates. Then, starting from two single BHs in KS coordinates (and different coordinate 
origin) a Lorentz boost is applied to each of them and then, by identifying the coordinate systems, a superposed 
metric is constructed. This metric is no longer an exact solution of the Einstein field equations (the constraint are 
not satisfied), but it is an approximation which improves as the BH separation increases. This type of data has been 
shown to work well in numerical relativity simulations (see, e.g. j24j). It has the advantage that one can consider Kerr 
BH initial data without the need to worry about spurious radiation that may originate in initial data constructions 
based on conformally-flat slicings. 

The evolution of the data is performed in analogy to that of BH-data, with one exception: the addition of matter 
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source terms on the right hand side of the evolution equations. We thus obtain a fully non-linear evolution of the 
geometry of the spacetime. In contrast, we substantially reduce the degrees of freedom of the matter sources. We will 
discuss in detail in this section, how this enables us to obtain an approximate evolution of the matter data. 



A. Description of the Spacetime Geometry 



We begin our description of the mathematical model with the spacetime geometry. In this work we follow the 
standard 3+1 splitting of Einstein's equations based on the Arnowitt-Deser-Misner (ADM) formulation [2^ (see [2|| 
llll for details). There it is assumed that the spacetime has topology Ix S, and hence can be foliated by the level 
surfaces E(£) of a function t(x^). The unit normal to the hypersurfaces £(i), — —aV^t, satisfies g lll ,n^n u = —1 . 
It can be used to decompose the four-metric according to 



(1) 



where W v is the projector orthogonal to (h} L v n v = 0). A vector field t^ threading the spacetime is introduced, 
fj,t = 1 , in terms of which we can introduce the lapse function and shift vector through the relations a — —t^n^ 



and = t» 



respectively. Then, using a system of spatial coordinates {x 1 } adapted to the hypersurfaces £(i), 



the metric takes the well-known form: 

ds 2 



-a 2 dt 2 + ^j{dx l + i3 i dt)(dx j + ftdt) 



(2) 



where 7^ = hij are the components of the spatial metric of the hypersurfaces The extrinsic curvature of E(f) 

is defined as 



K flu — 2 ^nhfiu ; 

where £ denotes the Lie derivative operator. We thus obtain 

where we have introduced the time derivative operator: 

d t = d t - £p . 



(3) 
(4) 
(5) 



Because we are dealing with a physical system that contains matter, we also need to apply the 3+1 splitting to the 
energy-momentum tensor T^ v . This can be done by decomposing T Miy with respect to the unit normal n^: 



T» v = En»n v + W + 2 ,]^n v) + IP" , 
where the different quantities that appear in this expression are defined as follows: 



E 



P — -h T^ v 

= -W v T VT n T , 



n 



flu 



— h h T ra — Ph 



Applying the 3+1 splitting to Einstein's equations we obtain the Evolution equations: 



dtlij = 
dtKij -- 

and the Constraint equations: 



-2aK i3 , 
-DiDjd - 



Ri 



KKi. 



2KikK k j 



n, 



1 



(E - P) 7 , 



H = R + K 2 - K l3 K i] - 2E = 0, 
V 1 = Dj{K l i - Y j K) - J 1 = . 



(6) 



(7) 
(8) 
(9) 
(10) 



(11) 
(12) 



(13) 
(14) 



Here, Di is the covariant derivative associated with the 3-metric 7^, iiy is the three-dimensional Ricci tensor, 
R = 7 y Rij the Ricci scalar, and K is the trace of Ky. Equations 1|11I12|) and (|13I14|1 constitute the basic 3+1 setup 



4 



and are known as the ADM equations [25j . We are going to use a modification of this formulation due to Shibata and 
Nakamura 0|, an d Baumgarte and Shapiro [ll|, now known as the BSSN 3+1 formulation of Einstein's equations, 
which has been found to result in vastly improved stability properties in numerical simulations compared with the 
ADM equations. The way in which the BSSN formulation is introduced from the ADM formulation is as follows: We 
introduce new variables based on a trace decomposition of the extrinsic curvature and a conformal rescaling of both 
the metric and the extrinsic curvature. The trace-free part Aij of the extrinsic curvature is defined by 



A; 



(15) 



Next, one introduces a conformal metric 7y in terms of the physical metric by 

7y = i>%j • 



(16) 



The value of the conformal factor ip can be fixed by demanding the determinant of the conformal metric 7^ to be 
unity: 



^ = 7 1/12 . 



lij = <F 4 7y = 7" 1/3 7« . 7=1, (17) 
where 7 and 7 are the determinants of 7y and 7^ respectively. Then, instead of 7^ and K^j we can use the variables 



= In ib = — In 7 . 

v 12 



K = 7yiT- 



Jij = 



7ij 



(18) 

(19) 
(20) 
(21) 



where %j has determinant 1 and has vanishing trace. Moreover, we also introduce the following conformal 
connection functions 



r 



7 



(22) 



where T l jk denotes the Christoffel symbols associated with the conformal metric. To obtain the second equality we 
have to use the fact that the determinant of the conformal spatial metric, 7, is unity. 

The variables used in the BSSN formulation are: (j), K, 7^, Aij, and P. The set of evolution equations for these 
variables that we are going to use, and which can be derived from the previous relations, is given by 



1 



d t jij = -2aA. 
d t <b 

B t K = -D l D t a 



— aK . 
6 



(DiDja 



TF 



(Rff - 8TrUij) + a(KAij - 2A ik A k j) , 



AiiA ij + -K 2 
y t 3 



4ir(E + 3P) 



d t r = f3 3 djr - r'djF - 2A ij djo 

+ j^djdkF + W j djd k f3 k + 



2a 



2~ 
3 



rdjft 



)(« 



o 



(23) 
(24) 

(25) 
(26) 

(27) 



djp 



Here the superscript TF denotes the trace- free part (with respect to the metric 7ij). In this context, it is important 
to note that <fr is defined in terms of the determinant of the 3-metric 7^, which is a scalar density of weight 2, and that 
7ij and A^ are tensor densities of weight —2/3 . Finally, x is a free parameter [it multiplies a quantity that vanishes 
at the analytic level by virtue of the definition of the quantities T l l|22|)] that has been set to 2/3 for all simulations 
reported here. This value has been suggested by numerical experiments. 



B. BH description: Kerr-Schild form of a Schwarzschild BH 



Although spinning BHs are likely to be more interesting from an astrophysical point of view, especially in the 
case of extreme-mass-ratio binaries, in this work we restrict our attention to non-rotating Schwarzschild BHs. We do 
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note, however, that for the case of KS coordinates, the fundamental properties of spinning and non-spinning black 
holes are the same from a numerical point of view (see, for example, the successful evolutions of single Kerr BHs 
obtained in [13 )• This feature represents a key motivation in our choice of coordinates and is not obviously satisfied 
by coordinates which are not adapted to the spinning cases (for instance, all the coordinate systems that assume an 
initially conformally-flat spatial geometry). 

In this section we describe the contribution to the initial data of the BH, that is we need to prescribe initial data 
for the functions ((f), K, 7^, Ay, V) . These arise from the KS form of the metric, which can be conveniently written 
as the sum of a flat space metric plus a term that factorizes as the product of a light-like vector with itself according 
to 

ds 2 = (7?^ + 2Ht fi £„)dx> 1 dx v . (28) 
Here r)^ u = diag(— 1, 1, 1, 1) is the Minkowski metric in Cartesian coordinates {t, x 1 } and £ M is a null vector field 

I dx" = -dt - —dx i , r 2 = SijxV , (29) 
r 

where Xi = x l . This light-like vector field has been chosen so that it corresponds to ingoing light rays. Finally, H is 
a scalar given by 

jj _ -Mbh ^ (30) 

r 

where M BH is the BH mass. 

The metric quantities are then independent of the time t which slices the spacetimc 

= Y+2H ' R" 1 = T+2H e ' (31) 
7 g H = % + 2H£ i £ j , 7 & = & Y^m ei3 ' (32) 
7 bh = det( 7 f/) = 1 + 2H = a-l . (33) 

C. NS matter description: The Tolman-Oppenheimer-Volkoff stellar model 

The main approximation in our description of the BH-NS binary system lies in the description of the NS. The 
main goal of this approximation is to provide a framework in which one avoids solving the hydrodynamical equations 
governing the NS matter fields while at the same time being able to have a reasonable approximation of the main 
dynamical aspects of the evolution of such a system. The implementation of this idea requires two ingredients: (i) 
The prescription of the profiles for the different matter variables and the construction of initial data for the geometry, 
(ii) The evolution of the matter profiles. Regarding the second point, in this paper we take the simplest description, 
namely, to consider rigid profiles that follow a prescribed trajectory. In subsection III El we briefly discuss how this 
description can be refined in future work to obtain physically more realistic models. 

With regard to the first point, the matter profiles for the NS matter variables are obtained by using an exact 
solution representing an isolated NS star, more precisely a TOV model. With regard to the eventual superposition 
with the BH, we will apply two modifications to the standard TOV metric: i) a transformation to KS-like coordinates, 
and ii) apply a boost transformation to account for the motion of the star relative to the black hole. Our starting 
point is the TOV metric in Schwarzschild coordinates 

g a/3 dx a dx p = -A 2 (R)df 2 + B 2 (R)dR 2 + R 2 dn 2 , (34) 
where dfl 2 = d9 2 + sin 2 9dip 2 . The energy momentum tensor of a TOV star is that of a perfect fluid and is given by 

T a p = (p + p)u a u p +pg a/3 , (35) 

where u a = (—A, 0, 0, 0) is the 4-velocity of the fluid elements. It is convenient to introduce the mass function m(R) 
by 

R- 2 = l-2m/R. (36) 
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Then, the TOV equations are: 



A' m + 4irp R 3 
1 ~ R(R~ 2m) 
m' — inpR 2 , 



(37) 
(38) 



where a prime denotes differentiation with respect to R . We are going to consider the case of a barotropic Equation 
of State (EoS), specifically a polytropic one: 



P = k P r , r = i + ~ 

n 



(39) 



where the polytropic coefficient k and the polytropic index n are constants. Equations I|37I38|I determine a one- 
parameter family of TOV solutions, which can be parametrized by the value of the central density p c = p{R = 0) . 
The sound speed of the matter distribution is given by c 2 = dp/ dp . 

Bearing in mind the eventual superposition of the BH and NS metrics in analogy to the construction of KS 
binary BH data, we are going to apply coordinate transformations to the TOV solution that resemble the coordinate 
transformations used in the BH case. In the case of BHs, the first coordinate transformation aims at changing the 
time coordinate so that it is adapted to light-like geodesies, typically ingoing ones. For a Schwarzschild BH, this 
corresponds to the KS coordinate system. In the case of a NS we follow the same idea. A straightforward calculation 
shows that after applying the coordinate transformation 



df = dT+{B/A-l)dR, dR = dR. 



(40) 



the ingoing radial light-like geodesies are given by T + R . This transformation, however, is not suitable for a 
three-dimensional numerical treatment because it leads to a coordinate singularity at the centre of the star when we 
transform to Cartesian coordinates [by using the transformation given below in equation l|44|) ]. It is possible, though, 
to find a coordinate transformation analogous to l|4(J|l that satisfies all regularity requirements at the stellar centre. 
For this purpose we define the new coordinate time T by 



4 r = *r + (*--L) <i fi. 

The resulting line element is then given by: 

-A 2 dT 2 + 2AB{1 - B- 2 )dTdR + (2 - B- 2 )dR 2 + R 2 dn 2 . 



ds 2 



(41) 



(42) 



By virtue of Birkhoff's theorem, the geometry exterior to the star is necessarily described by the Schwarzschild metric. 
In order to combine the interior and exterior metric, we need to match them in a smooth way at the stellar surface. It 
is a remarkable property of the interior metric (|42l) that it still smoothly matches to the exterior (KS) solution and thus 
leads to coordinates adapted to the ingoing null geodesic structure outside the star. This follows straightforwardly from 
performing the substitution A 2 = B~ 2 = 1 — 2M NS /R, where M NS is the NS mass. In conse quen ce, the interior and 
exterior geometries match smoothly (in the sense of the Lichnerowicz junction conditions ,29), 30]), at a hypersurface 
R = i? NS = constant , provided the TOV model satisfies the following conditions 



A 2 \r=r nr = 1 - 2Af NS /i? NS , p(i? NS ) = 0, (43) 

r=h n „ = 1 - 2M NS /i? NS (or equivalently AB\ R=R = 



-21 



which have the obvious consequences: m(R NS ) = M NS and B 
1), and clearly i? NS represents the NS surface. 

Next, we apply the coordinate change from spherical (T,R,9,ip) to the Cartesian-like coordinates (T, X 1 ) . These 
coordinates are related by 



(X 1 ) = (X, Y, Z) = {R sin9 cos tp, R sin 9 simp, Rcost 



The resulting line element is 



ds z = —A dT + 2AB(1 - B- z )dT—dX 

R 



(1-B 



-2\XiXj 



R 2 



dX l dX 3 



(44) 



(45) 



where R is considered a function of X 1 given by R 2 = 5ijX l X^ , and Xi — X 1 . At this point, we can see that B — > 1 
as R — > , i. e. the metric is regular at the origin X 1 = of the stellar object. This would not have been the case for 
the coordinate transformation l|40|l which would have lead to a cusp at the origin of the star. 
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The ADM metric variables associated with the slicing {T = constant} are: 

A 2 B 2 



Pc = AB 



2-B- 2 ' 

1 - B- 2 X 1 



B- 2 R 



(46) 
(47) 
(48) 



where we have used the tag V to distinguish the quantities associated with this particular slicing. We now address 
the next step in the construction of the superposed geometry, the Lorentz boost. We first note, that the NS metric 
(|45|l cannot be cast in the KS form (|28fl. which is preserved under Lorentz transformations. While the NS metric can 
still be written in the form t]^ + the resulting H^ v cannot be factorized in terms of a null vector field, where 
Hn V ~^ as R — > oo. On the other hand, the form rj^ u + H^ v is still Lorentz invariant. We thus apply a coordinate 
change from the Cartesian-like coordinates (T, X 1 ) of Eq. i|45[) to a different set of Cartesian-like coordinates (t,x l ) 
via a Lorentz boost. These two sets of coordinates are related by the following expressions: 



T 

X 1 



1 [t-8 lJ v\x 3 -z 3 )} 



-jtv 1 



(x 3 — z 3 ) , 



(49) 
(50) 



where v l are the components of the boost velocity (vi = v l ), v 2 is 8ijV l v 3 , 7 is (1 — v 2 )^ 1 / 2 , and z l are the spatial 
components of the compact object's location. In the course of the evolution, the position vector z % needs to be evolved 
by solving the equations governing the external dynamics of the NS. In this paper it is going to be a prescribed 
trajectory. 

After carrying out the transformation (|49I50|I the ADM metric variables associated with the slicing {t = constant} 
can be written as 



i 2 



Q 

B l = - 



a 2 cl ~ 2 [(1 - vrfl) 2 



7 



7.. 



1 + 7 



7 a NS 



"o7c v i v j] 

: 2 {l-vifc) ( -K 

7 



1 + 7 



-v 



1 + 7 



'7 



7w A' 



+ 7o«j 
7 



1+7 



ft 



1 -v l 



1 + 7 



VjV 



i u 3 1 



(51) 
(52) 

(53) 



where a a , j3 l a , and 7° are given in equations l|46|l - H48(l in terms of the coordinates (T, X 1 ) . With expressions 151ll - 153() 
we complete the geometric description of the NS. They will be used below to construct the initial data by superposing 
the BH and the NS solutions. 

We denote the matter source terms associated with the foliation {t = const.}, by (E NS , P NS , J^ s , n^ 3 ). These can 
be obtained from the expressions (f7 | - (|lt)|l . bearing in mind that the normal to this foliation is n^dx^ = —a NS dt and 
the components of the fluid four-velocity are 



where u T = A 1 . Then, the result is 



(54) 



j' 



n NS 



ip+p){— 



(P + P)^ hTiJt - hff-rki) (A- + v k )(p l NS + v l ) . 



(55) 
(56) 

(57) 
(58) 



where a NS , /9Jj S , and 7^ s are given by l|51 [l -l|5H jl . and p and p are the TOV energy density and isotropic pressure. 
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D. Superposition of the BH and NS: Initial Data and BSSN Source Terms 



We now have all the necessary ingredients to superpose the BH and NS spacetimes given above. Because of the 
presence of matter terms, this process involves two steps: (i) the construction of the geometric quantities associated 
with the BH-NS system and, (ii) the construction of the matter source terms associated with that geometry. 

First, we superpose the geometries in complete analogy with the BH binary case (cf. Let us consider, for 

this purpose, the BH metric given by expressions (|28|) - (|3U[I and the NS metric given by expressions (|51fl - l|53l) . By 
identifying the two spacetimes' coordinate systems, the spatial metric and extrinsic curvature of the superposition are 
given by 



lij 



™ BH _|_ ^ N .s 

K 



(59) 
(60) 



Moreover, the lapse and shift associated with the superposition are given by the following expressions: 

a = (a^ + a-' - l)" 1 / 2 , (61) 
p i = y^BH + ^s) _ (62) 

Second, we need to obtain from the NS description the expressions for the matter sources that appear on the right- 
hand sides of the BSSN equations 1(25 ^ -^25 )1 . To that end, we consider the following form of the energy-momentum 
tensor 



= (p + pWW+pgr , (63) 

where g^ v is the inverse of the global metric given in equations l|59l) - (|62|) . p and p are the TOV energy density and 
isotropic pressure written in the coordinates (t, x l ) of equations (|49I50|I . and C/ M = U t {l,v % ) is the fluid velocity. The 
time component U i of the fluid velocity is determined in terms of the boost velocity v l and the global geometry by 
imposing the normalization condition g a f3U a U 13 = — 1. The resulting expression is 



m 2 



a 2 



lkl (P k +v k )(p l +v l ). (64) 



Following the same procedure as above, we can compute the matter quantities (E,P, J l ,Hij) associated with the NS 
in the superposed BH-NS spacetime. We note that the normal to the foliation is now n^dx^ — —adt. We thus obtain: 

E = (p + p)(aU t ) 2 -p, (65) 
P = p+±( p + p)(U t ) 2 7ij (p i + v i )((3i +vj)=p-±(E-p), (66) 



J 1 = (p+p)a(U t ) 2 (^+v t ), (67) 
% = (p+p){U t ) 2 ( 1 . lkl]l -\ ll3lkl ){p k +v k ){p l +v l ) 1 (68) 

where a , f3 l , and 7^ are given by l|61[) . (|62() . and (|59() respectively, and p and p are the TOV energy density and 
isotropic pressure. These expressions are very close in form to the expressions we obtained for a single NS, the main 
difference being that the geometric objects now refer to the superposed spacetime instead of the TOV spacetime. 

It is clear that these initial data only satisfy the constraints in the limit of infinite BH-NS separation. The constraint 
violations inherent to this construction, however, are expected to be small. If, for example, one considers this type of 
superposed data as the b ackg round metric fields for the conformal transverse-traceless method of Lichnerowicz, York 
and others EU- l3fij | - the conformal factor obtained from solving the Hamiltonian constraint is close to unity 

throughout the whole spacetime 23] . Furthermore the quality of the data can always be enhanced by increasing the 
distance between the BHs. Finally, this data is believed to avoid spurious radiation contamination, and in contrast 
to conformally-flat data, it is possible to introduce Kerr BH solutions in a straightforward way. 

Bonning et al (2^ have further shown that the superposed binary BH initial data lead to the correct (Newtonian) 
binding energy in the Newtonian limit. We can ask the same question about the prescription we have presented 
for BH-NS binaries. The answer is negative, and, in essence, the reason is that the way in which the NS has been 
introduced does not account for the deformations caused by the BH gravitational field. In the case of BH binaries 
the binding energy is defined as E b = £^ - £^ H - _Eg H , where the total ADM energy is Mg H + M| H and the 
individual energies are defined in terms of the apparent horizon masses: E^ H — Mf H (I = 1,2). In the Newtonian 
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limit (GM/(c 2 £) <C 1 , where i denotes the coordinate separation) they are given by [23[: E\ n « Mj* H [l + M| H /(2^)] 
and £| H s» M| H [1 + Mg H /(2£)]. We thus get the familiar expression £5 « -Mg H Mj H /£. If one of the objects is a 
NS, all that changes, is the individual energy E NS associated with the star. It is now obtained from the integration of 
the NS energy density over the initial slice. Including the contribution of the BH to the volume element, the result 
is E NS ss M NS (1 + M BH /£) — E^", where E^l 1 denotes the NS self-binding energy. Therefore, we do not recover the 
Newtonian expression. The reason for this lies in the way in which the NS has been introduced in the superposition 
and the fact that the definition for its individual energy, E ns , does not include the gravitational influence of the BH. 
It turns out that this way of superposing the NS is closer to the case of a test body. Indeed, if we look at the binding 
energy by considering the NS as a test body, we find that Eb — M NS (1 — E NS /M NS ), i. e. the usual definition for test 
masses, and we obtain the correct Newtonian expression. 

We could modify our description by performing an additional coordinate change to the NS model, prior to the 
superposition. This would result in extra parameters, which may be fixed so that the deformations of the NS produced 
by this coordinate change will account for the effects caused by the presence of the BH. While a detailed study of 
such modifications is beyond the scope of this paper, it will be interesting to study their effect in future work. 

E. On the Equations of Motion for the NS 

In a completely general setup the equations of motion for the matter fields associated with the NS just follow from the 
local conservation of energy and momentum, i.e. V^T M1/ = , which leads to the equations of magneto- hydrodynamics. 
In this paper, however, we are advocating an approximate approach that avoids solving hydrodynamical equations 
by reducing the number of physical degrees of freedom to a finite number, the hydro without hydro approach. More 
precisely, the idea is, that we can approximate the matter description by a simplified model that involves a finite 
number of parameters, {A 7 }, so that the associated energy-momentum distribution is determined in terms only of 
these parameters and the trajectory of a certain center of mass, z^, that is, = T^ v \z p ; A 7 ] . Then, in order to 
update the matter sources that enter in our BSSN evolution equations we merely need to evolve the parameters {A 7 } 
(referred to in this work as the internal motion) and the components of the trajectory z M (the external motion). In 
practice, evolution of these parameters is obtained from ordinary differential equations. This is a direct consequence 
of going from an infinite number degrees of freedom, as described by partial differential equations, to a finite number, 
described by ordinary differential equations. 

For the case of extreme-mass-ratio binaries, the structure of the small object may be neglected as a first approxi- 
mation, so that we do not need to consider the interior motion. The external motion, in turn, is given by the solution 
of the geodesic equations in the numerically constructed spacetime. In order to obtain a more realistic description of 
BH-NS binaries with similar masses, it is important to take into account the internal motion or, in other words, we 
need to allow for deformations of the NS associated with tidal forces arising from the presence of the BH. To that 
end, an interesting approach is to consider the so-called affine stellar model introduced by Lattimer and Schramm j3?J 
and further developed by Carter and Luminet [3^, |3!| and generalized for curved spacetimes in [40|. As it stands 
now, however, the affine model is applicable only to slowly varying general relativistic spacetimes |4Ctl4lll4^ |. where a 
stationary background can be identified. This is an assumption that cannot be adopted in our framework, since we are 
interested in situations that involve strong and dynamical gravitational fields. Therefore, the affine model should be 
generalized appropriately for dynamical spacetimes. A particularly interesting way of doing this is based on Dixon's 
work [43L fiiL I45ll4q | (see also H^)- The idea is to use the multipole moments associated with the energy-momentum 
tensor to construct equations of motion for the center of mass, i.e. 2^, and for the extra parameters that describe 
the NS internal motion. In any case, the affine model has been used successfully to study the tidal interaction and 
disruption of stars by supermassive BHs isolated rotating stars jl^, stars in binary systems jl^ and, more 

recently, gravitational signals arising from tidal interactions [50l |. Another semi-analytic approach is the ellipsoidal 
energy variational or Roche- Riemann model |5lj and its relativistic generalization |52j which is formally equivalent to 
the hydrostatic limit of the affine model. 

In this paper, our goal is to investigate the stability of the numerical relativistic evolutions under the presence 
of moving matter distributions. To that end, we adopt the simplest possible evolution for the NS, namely a fixed 
trajectory around the BH. 

III. COMPUTATIONAL FRAMEWORK AND NUMERICAL COMPUTATIONS 

The numerical results presented in this work have been obtained with the 3D numerical relativity Maya code which 
has been described in detail in [l3L ITij . The code uses the Cactus computational toolkit for parallelization, 
data input/output and horizon finding. Inside the Cactus environment mesh refinement is provided by the Carpet 
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TABLE I: Parameters and resulting mass and density for the six stellar models discussed in this work. 



Model 


k 


n 


-Rns [Mbh] 


M NS [M bh ] 




pc [M B i] 


1 


0.831 


1 


1 


lO" 1 


IO" 1 


7.9 10" 2 


2 


0.652 


1 


1 


10~ 2 


10~ 2 


7.9 10~ 3 


3 


0.638 


1 


1 


10 -3 


10 -3 


7.9 10~ 4 


4 


0.00831 


1 


0.1 


io- 2 


io- 1 


7.9 


5 


0.00652 


1 


0.1 


10 -3 


10 -2 


7.9 10" 1 


6 


0.00638 


1 


0.1 


1Q -4 


10" 3 


7.9 10~ 2 



package [T^. With regard to earlier versions of the Maya code used for BH simulations, the evolution of extreme- 
mass-ratio binaries in the framework of the geodesic approximation method has made necessary some additions and 
modifications to the code. Foremost, these are the implementation of dynamic (moving) mesh-refinement and the 
addition of the matter source terms on the right hand sides of the Einstein field equations. Both will be described in 
more detail as we discuss the numerical simulations performed in this work. 



A. Description of the Numerical Implementation 

We first address the calculation of the TOV profiles. The TOV equations lpT7 jl -l(55 )l constitute a boundary value 
problem for the variables A(r), m(r) and p(r). This system of equations is closed by an equation of state p = p(p), 
chosen to be of polytropic type in this work. The boundary conditions are given by m = at the center as well as 
p = and A = y/l — 2M NS /i? NS at the stellar surface. We solve the resulting system of ordinary differential equations 
using a standard relaxation scheme j^. In practice, we calculate the TOV profiles during the initialization phase of 
the code and store them as functions of the radial variable r. The values required on the three-dimensional Cartesian 
grid are then obtained from third order polynomial interpolation. 

The discretization of the BSSN evolution equations (|23I28(I has been implemented using centered second order 
stencils for spatial derivatives, except for the advection terms of type f3 k dk, for which second order accurate upwinding 
operators have been used. The integration in time is done using an iterative Crank-Nicholson scheme. As in simulations 
of vacuum spacetimes we have treated the spacetime singularities associated with the BHs by excising an area of finite 
size inside the apparent horizons (see |l3L ll 1| for the details of the implementation of this technique) . 



B. Analysis of the Initial Data 



In order to study the properties of the initial data proposed in this paper we have considered six different TOV 
models for the NS. The physical features of these models (the polytropic coefficient fc, the polytropic index n, stellar 
radius i? NS and mass M NS , compactness ratio M NS /i? NS , and central density p c ) are shown in TableQ] Notice, that the 
compactness ratios are small and hence some of the models would describe stars much less compact than astrophysical 
NSs. 

From the values in Table |U we find that the central density does not only increase with the compactness of the 
stellar model, but also when we make the star smaller, keeping constant the compactness ratio (see, for instance, 
models 1 and 4). This is an expected result: if we reduce the mass M NS while keeping the compactness ratio M NS /i? NS 
constant, the density behaves like 

M ns J_^L _ const 
P ~ 7? 3 ~~ M 2 7? 3 ~~ M 2 ' 1 ' 

NS J "NS NS 1U NS 

That is, the density scales as the inverse of the square of the mass. The values in Table H] confirm this relation. 
From the computational point of view, this feature represents a considerable challenge. For more extreme mass ratios 
(smaller M NS ) we do not only need to decrease the grid spacing linearly with the size of the star. Further increase in 
resolution is necessary to adequately resolve the steep gradients resulting from the increasingly larger matter densities 
and curvature encountered for smaller values of M ns . Even with the availability of mesh refinement, as discussed 
below, the computational demands for simulating low mass compact sources quickly become prohibitively costly. It 
is for this reason that we also consider the less compact models listed in Table [i] We emphasize, however, that all 
models under consideration here are orders of magnitude more compact than those used in similar studies in the 
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TABLE II: Structure of the refinement levels as used in the study of the initial data. The two values for the grid spacing 
correspond to the coarse and fine resolution in the convergence analysis. The boundaries in the z-direction are z m i n = and 
Zmax = y m ax on all levels and have not been listed explicitly. Levels 7-9 are only used for stellar models 4-6. 



Refinement 
Level 



hi [Mbh] h 2 [M BH ] 



Grids centered around the BH 

[M BH ] ■''max [M BH ] 1/min [M BH ] V max [Mb H I 



Grids centered around the NS 
ain [M BH ] Zmax [M BH ] J/min [M BH ] V [Mb Hi 



2.88 
1.44 
0.72 
0.36 
0.18 
0.09 
0.045 
0.0225 
0.01125 



1.92 
0.96 
0.48 
0.24 
0.12 
0.06 
0.03 
0.015 
0.0075 



-103.68 
-51.84 
-25.92 
-12.96 
-4.32 



103.68 
51.84 
25.92 
12.96 
4.32 



-103.68 
-51.84 
-25.92 
-12.96 
-4.32 



103.68 
51.84 
25.92 
12.96 
4.32 



6.48 


10.80 


-2.16 


2.16 


7.56 


9.72 


-1.08 


1.08 


8.1 


9.18 


-0.54 


0.54 


8.37 


8.91 


0.27 


0.27 


8.505 


8.775 


0.135 


0.135 



literature [2 11) . Below we will discuss the numerical demands arising from high compactness more quantitatively, 
when we discuss the constraint violations of the superposed data. First, we test the numerical implementation of the 
boosted TOV solution by studying the convergence of the constraints in the absence of a BH. 

The initial data constructed according to the procedure described above results in exact solutions of Einstein's 
equations in various limits. In particular, by setting the mass of the BH to zero we recover the solution of a boosted 
TOV star in KS-like coordinates. We emphasize, that, in the absence of the BH, there is no fundamental length scale 
M BH , so that the numerical values obtained for the NS no longer have the unambiguous meaning as before. A key 
motivation for this study, however, is to probe the suitability of the numerical framework used for the simulations. For 
this purpose we construct the numerical grid using the same numerical values as if a BH of mass M BH were present. In 
this work we limit our discussion to non-spinning BHs and stellar models. In consequence, the configurations under 
study are inherently symmetric about the orbital plane. It is sufficient, for this reason, to evolve data only in the 
bitant z > and impose symmetry conditions in the orbital x, y plane. 

In order to achieve maximum resolution at moderate computational cost, we use Fixed Mesh Refinement for all 
simulations in this work. Inside the MAYA code, Mesh Refinement is implemented via the CARPET package 0,^3- 
In contrast to the static fixed mesh refinement used for the head-on collisions in [24| the type of scenario investigated 
here will eventually require dynamic mesh refinement to accommodate the orbital motion of the stellar object around 
the BH. We return to this issue later, when we will discuss time evolutions. For the initial data under consideration 
here, we restrict our attention to the initial setup of the refinement levels. The initial position of the star in this case 
is x — 8.64, y = z = and it has an initial boost velocity given by v x = v z = , v y = 0.34. For the stellar models 
1-3 of Table H] we use a set of 5 nested boxes around the BH plus a set of two nested refinement components centered 
around the NS. For models 4-6 we add three further refinement levels centered on the NS. The exact specifications of 
the refinement levels are given in Table ITT1 

For the convergence analysis we have calculated the data for two different grid spacings hi and li2 — h\j '1.5 as listed 
in the table. Because an isolated boosted TOV star is a solution of the Einstein equations we expect the Hamiltonian 
and momentum constraints, given by equations I|13ll4f) . to converge to zero at second order. Numerically this implies 
that the constraint violations obtained for the high resolution calculation should be 1.5 2 smaller than those obtained 
at coarse resolution. 

Figure ^ demonstrates second order convergence for the Hamiltonian and the norm ||.M|| = yjjj Mj of the 
momentum constraints along the x-axis, obtained for the most compact model of Table [IJ model 4. We find similar 
results for all constraint functions on arbitrary axes. A similar analysis at coarser resolution (1.5 hi instead of /12), 
however, gave less satisfactory results, indicating that the resolutions used here are necessary for model 4. For the 
other, less extreme, models we find correspondingly less strict requirements on the resolution to achieve second order 
convergence. 

We now turn our attention to the binary data obtained for non-vanishing mass of the BH. In the continuum limit the 
resulting data are a solution of the Einstein equations only as the separation of the two objects becomes infinite. We 
therefore do not expect the constraint violations to converge to zero. On the other hand, the constraints will never be 
satisfied exactly in a numerical simulation because of the finite numerical accuracy. In order to assess the significance 
of the constraint violations inherent to the superposed data relative to the limitations in numerical accuracy we follow 
the same approach as in |l4l ] and study the constraint violations at different resolutions. If the constraint violations 
decrease at about second order for higher resolutions we conclude the inherent constraint violations (due to the fact 
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FIG. 1: Constraint violation of the Hamiltonian and the norm of the momentum constraint along the 2>axis for stellar model 
4 in the absence of a BH. The constraints obtained at high resolution have been amplified by 1.5 2 . 



that we are not solving them) to be insignificant relative to the purely numerical error. 

For this purpose we show in Figure |3 the convergence of the Hamiltonian constraint and the norm ||A4|| of the 
momentum constraints along the x-axis. For clarity we have used a logarithmic scale for the distance from the origin 
in both directions. As in the case of an isolated star, we find that the constraints converge away at second order 
for the resolutions used here. The only exception is the region of the star, where the non-linear interaction between 
the gravitational fields of the two objects is strongest and we expect the constraint violations inherent to this data 
construction to be largest. It is these constraint violations which represent the approximative nature of our approach 
to study this type of binary configurations. 

Figure|21also illustrates the high computational demands for simulations of more compact objects. The two models 
(3 and 4 of Tabled used in the figure represent the least and most compact stars. Clearly the constraint violations 
found for model 4 in the upper two panels are very large compared to the average violations near the BH. This 
is an artifact of the high densities and curvature encountered near the small star in combination with resolutions 
currently available with our computers. In contrast the constraint violations found for model 3 in the lower panels is 
substantially smaller than those observed near the BH. We emphasize that this issue is separate from the convergence 
properties mentioned above. In fact, the convergence properties found for model 4 are closer to second order than 
those found for model 3. We attribute this to the fact that the constraint violations near the very compact object 4 
are dominated by the discretization errors arising from the insufficient resolution near the stellar center. For these 
reasons we study in the next subsection evolutions of the less compact model 3 only, for which we are confident to 
achieve sufficient resolution, even inside the star. 



C. Time evolutions 



In this section we present time evolutions of the initial data presented above. Following our previous discussions, 
we focus on the implementation of the different elements of this type of simulations in the framework of moving 
refinement components with emphasis on the stability of the evolutions. 
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FIG. 2: Constraint violation of the Hamiltonian and momentum constraints along the :r-axis for stellar model 4 (left panels) 
and 3 (right panels). The constraints obtained at high resolution have been amplified by 1.5 2 . 



TABLE III: Structure of the refinement levels as used for the time evolution. 



Refinement Level 


h [M BH ] 


Imin [M BH ] 




J/min [M BH ] 


J/max [M B h] 


I 


4.00 


-144.00 


144.00 


-144.00 


144.00 


2 


2.00 


-72.00 


72.00 


-72.00 


72.00 


3 


1.00 


-36.00 


36.00 


-36.00 


36.00 


4 


0.50 


-18.00 


18.00 


-18.00 


18.00 


5 


0.25 


-12.00 


12.00 


-12.00 


12.00 


6 


0.125 


5.00 


9.00 


-2.00 


2.00 



To this end, we have evolved the superposed data of a BH of mass A/ BH = 1 and stellar model 3 as listed in Tabled 
Because of memory restrictions of the available computational resources we have evolved this binary system using a 
distance d = 7 M and a grid setup slightly different from that of Tabled The grid specifications used in the evolution 
are listed in Table ITTT1 Of the refinement levels listed there, number 1 to 5 remain fixed throughout the evolution. 
Level 6, however, is required to follow the motion of the stellar object and thus needs to move. For this purpose we 
use the regridding operation implemented inside Carpet. This operation allows for interpolation of function values 
from coarser onto finer refinement levels. This operation is regularly performed in mesh refinement simulations of 
Berger-Oliger type to provide boundary conditions for the inner levels, the so-called prolongation. In case of a moving 
high resolution grid component, the same operation is used to fill the new points on that component with valid data 
from the coarser levels. Inside the Maya code we use this feature by adjusting the specifications, i.e. x m i n , x max , t/ m i n 
and ?/max, of refinement level 6 in the course of the evolution. 

We have implemented two alternative ways of controlling these values. The first uses the tracking of the BH 
singularity, as provided, for example, by an apparent horizon finder, and moves the center of the refinement component 
by a corresponding amount. Because this method is inherently restricted to BHs, we also allow for the center of the 
component to follow a user specified trajectory. This second method is realized for the case of the stellar object 
discussed in this work. The trajectory is the same as that used for the motion of the stellar object itself. In Fig. we 
show four snapshots of the ensuing evolution. The figure shows the trace of the extrinsic curvature obtained at times 
t = 0, 30 M, 70 M and 120 M. While the BH is represented by the large central throat, the stellar object manifests 
itself in the form of the small perturbation initially visible to the right of the BH (upper left panel). 

The simulations we have carried out last at least for an evolution time of 167 M, corresponding to about one and a 
half orbits, without showing signs of instability. Furthermore, no elaborate fine-tuning of the parameters was necessary 
to achieve evolutions. In conclusion, simulations on orbital timescales appear to be achievable rather straightforwardly 
using the approach discussed in this work. 
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FIG. 3: Snapshots of the time evolution of a NS model type model 3 orbiting a massive BH with orbital frequency lo = 
0.054 M™ 1 . 



IV. CONCLUSIONS 

In this paper we have introduced a framework to study BH-NS binaries using a fully non-linear evolution of the 
spacetime geometry and an approximate evolution of the NS matter sources based on the reduction of the degrees 
of freedom of the NS, so that hydrodynamical calculations are avoided. The remaining degrees of freedom, including 
the motion of the NS are governed by ODEs. We expect this framework to be useful to investigate certain dynamical 
regimes of the BH-NS binary system relevant for the point of view of gravitational-wave astronomy. 

The construction of the initial data is performed by generalizing the superposition technique of [2(j , as applied 
to BH binaries in KS form, to BH-NS binaries in KS-like coordinates. In order to make this approach functional 
in numerical relativity using three-dimensional Cartesian coordinates, we found it necessary to modify the standard 
coordinate transformation to KS coordinates inside the NS. In particular, we used a coordinate transformation that 
avoids the appearance of a coordinate singularity at the stellar centre, while preserving the light-like geodesic structure 
of the coordinates outside the star. 

We have studied the properties of the initial data thus obtained for a set of stellar models with different values of the 
compactness and mass ratios. An adequate description of extreme-mass-ratio binaries involving compact stars is only 
made possible by the use of advanced mesh-refinement techniques. We have demonstrated second-order convergence 
of the Hamiltonian and momentum constraints in the case of an isolated NS (corresponding to an infinite BH-NS 
separation) . For finite separations we still observe second-order convergence of the constraints over large parts of the 
computational domain. We thus conclude that these initial data sets are as suitable for the description of BH-NS 
binaries as the pure BH-BH analog previously discussed in the literature. 

Finally, we have investigated the feasibility of evolving such scenarios in time using modern methods of numerical 
relativity. In the numerical simulations considered in this work for testing purposes, we have restricted ourselves to the 
simple case in which the NS follows a prescribed trajectory. Foremost, our interest was focused on the stability of the 
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resulting simulations and the demands on computational resources. We have found that for sufficiently extreme mass 
ratios, long-term stable simulations on orbital time scales can be achieved straightforwardly using standard evolution 
schemes for solving the Einstein field equations. As the equal-mass limit is approached, the stability properties 
deteriorate because of the strong dynamics in the vicinity of the BH. In this limit, our approach no longer provides 
an adequate description of the matter sources, however. The scenarios of immediate interest for our approximation 
technique are thus within the range of capability of our evolution code. 

With regard to the compactness of the NS we find that the combination of extreme mass ratios with high com- 
pactness of the NS result in very steep gradients. These, in turn, demand a very high resolution and a large number 
of refinement levels. Currently we do not have the computational resources available to provide such resolution in 
long-term evolutions. We have therefore studied the time evolutions of less compact objects, but still orders of mag- 
nitude more compact than what has hitherto been considered using this type of approach. Such studies are well 
within the scope of current resources and can be refined in the future, by adding radiation-reaction effects to the 
dynamics and/or internal degrees of freedom of the NS to study different aspects of these systems in the context of 
gravitational-wave astronomy. 
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